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We treat the probability distributions for quadratic quantum fields, averaged with a Lorentzian 
test function, in four-dimensional Minkowski vacuum. These distributions share some properties 
with previous results in two-dimensional spacetime. Specifically, there is a lower bound at a finite 
negative value, but no upper bound. Thus arbitrarily large positive energy density fluctuations are 
possible. We are not able to give closed form expressions for the probability distribution, but rather 
use calculations of a finite number of moments to estimate the lower bounds, the asymptotic forms 
for large positive argument, and possible fits to the intermediate region. The first 65 moments are 
used for these purposes. All of our results are subject to the caveat that these distributions are not 
uniquely determined by the moments. However, we also give bounds on the cumulative distribution 
function that are valid for any distribution fitting these moments. We apply the asymptotic form of 
the electromagnetic energy density distribution to estimate the nucleation rates of black holes and 
of Boltzmann brains. 
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I. INTRODUCTION 



There has been extensive work in recent decades on the definition and use of the expectation value of a quantum 
stress tensor operator. When this expectation value is used as the source in the Einstein equations, the resulting 
semiclassical theory gives an approximate description of the effects of quantum matter fields upon the gravitational 
field. This theory gives, for example, a plausible description of the backreaction of Hawking radiation on black hole 
spacetimes 

However, the semiclassical theory does not describe the effects of quantum fiuctuations of the stress tensor around 
its expectation value. Quantum stress tensor fluctuations and the resulting passive fiuctuations of gravity have been 
the subject of several papers in recent years pm4, . However, most of these papers deal with effects described by the 
correlation function of a pair of stress tensor operators, and ignore higher-order correlation functions. 

One way to include these higher-order correlations is through the probability distribution of a smeared stress tensor 
operator. This distribution was given recently for Gaussian averaged conformal stress tensors in two-dimensional fiat 
spacetime [TS]. This result will be discussed further in Sec. IIB A recent attempt to define probability distributions 
for quantum stress tensors in four dimensions was made by Duplancic, et al [16) . However, these authors attempt to 
define distributions for stress tensor operators at a single spacetime point. Because such operators do not have well- 
defined moments, the resulting probability distribution is not well-defined. In our view, only temporal or spacetime 
averages of quantum stress tensors have meaningful probability distributions in four dimensions. Furthermore, these 
averages should be normal ordered, resulting in a zero mean for the vacuum probability distribution and a nonzero 
probability of finding negative values. None of these conditions are satisfied by the distribution proposed in Ref. |16]. 

The purpose of the present paper is to obtain information about the form of the probability distribution for averaged 
stress tensors in four-dimensional spacetime from calculations of a finite set of moments. This method was used in 
Ref. '15' to infer the distribution for ip^ , with Lorentzian averaging, where Lp is a, massless scalar field four-dimensional 
Minkowski spacetime. The result matches a shifted Gamma distribution to extremely high numerical accuracy. 



'Electronic address: "chris.fewster@york. ac . uk 
t Electronic address: f ord@cosmos.phy.tufts.edu| 
t Electronic address: romanSSccsu.edu 



2 



Unfortunately, the probability distribution of the smeared energy density for massless scalar and electromagnetic 
fields cannot be found so precisely. However, under certain assumptions to be detailed later, we are able to give 
approximate lower bounds and asymptotic tails for these cases, and to give a rough fit to the intermediate part of the 
distribution. 

An important point arises here. Throughout this paper, all quadratic operators are understood to be normal-ordered 
with respect to the Minkowski vacuum state. However, the smeared normal ordered operators are defined, in the 
first instance, only as symmetric operators on a dense domain in the Hilbert space (assuming a real-valued smearing 
function) and it is possible that there is more than one way of extending them to provide self-adjoint operators 
The operators of greatest interest to us are bounded from below on account of quantum inequalities (see Sect. II Al and 



so there is a distinguished Friedrichs extension (see Ref. [17], Theorem X.23), whose lower bound coincides with the 
sharpest possible quantum inequality bound. It is this operator that we have in mind when we discuss the probability 
distribution of individual measurements of the smeared operator in the vacuum state. The question of whether there is 
more than one self-adjoint extension, i.e., whether the normal ordered expressions fail to be essentially self-adjoint, is 
nontrivial and not fully resolved. Recent results (not, however, immediately applicable to our situation) and references 
may be found in Ref. [IF. If there are distinct self-adjoint extensions, their corresponding probability distributions 
will all share the same moments in the vacuum state. 

This links to the wider issue of whether or not the moments of the probability distribution determine the distribution 
uniquely. There is a rich theory concerning this question, which is reviewed in Ref. }19j. As will be discussed below, 
some of the moments we study grow too fast to be covered by well-known sufficient criteria (due to Hamburger and 
Stieltjes) for uniqueness. This does not prove that the distribution is nonunique (nor would the existence of distinct 
self-adjoint extensions prove nonuniqueness) and we have not been able to resolve the question of uniqueness. However, 
in Sect. |Vl| we prove that any probability distribution with the moments we find has a cumulative distribution function 
close to that corresponding to the fitted asymptotic tail. As various applications (see Ref. [20] and Sect. VII I depend 
only on the rough form of the tail, the possible lack of uniqueness is not as crucial as might be thought. Further 



discussion of this point can be found in Sect. VIII A 



II. REVIEW OF SOME PREVIOUS RESULTS 



Here we will briefly summarize selected aspects of two topics, quantum inequality bounds on expectation values, 
and known results for probability distributions. Both of these related topics are important for the present paper. 



A. Quantum Inequalities 

Quantum inequalities are lower bounds on the smeared expectation values of quantum stress tensor components 
in arbitrary quantum states |21H27) . In two-dimensional spacetime, the sampling may be over either space, time, or 
both. In four dimensions, there must be a sampling either over time alone, or over both space and time, as there are 
no lower bounds on purely spatially sampled operators Here we will be concerned with sampling in time alone, 
in which case a quantum inequality takes the form 

/OO 
f{t){T{t,0))dt>—^, (1) 
-OO 

where T is a normal-ordered quadratic operator, which is classically non-negative, and f{t) is a sampling function 
with characteristic width t. Here C is a numerical constant, typically small compared to unity, and d is the number 
of spacetime dimensions. 

Although quantum field theory allows negative expectation values of the energy density, quantum inequalities 
place strong constraints on the effects of this negative energy for violating the second law of thermodynamics [21) . 
maintaining traversable wormholes [35] or warpdrive spacetimes [301 ■ The implication of Eq. ([T]) is that there is an 
inverse power relation between the magnitude and duration of negative energy density. 

For a massless scalar field in two-dimensional spacetime, Flanagan j25j has found a formula for the constant C for 
a given f(t) which makes Eq. ([T]) an optimal inequality, and has constructed the quantum state in which the bound 
is saturated. This formula is 

where f{t) = T~^g{u) and u = t/r. This is the c = 1 special case of a general result for unitary, positive energy, 
conformal field theories in two dimensions, where c is the central charge, in which the left-hand side of (l2|) is multiplied 
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by c [55]. In four-dimensional spacetime, Fewster and Eveson P7| have derived an analogous formula for C, but in 
this case the bound is not necessarily optimal. 



B. Shifted Gamma Distributions 



Here we briefly recall the main results of Ref. [TS]. First, we determined the probability distribution for individual 
measurements, in the vacuum state, of the Gaussian sampled energy density 

Ttt{t,0)e-'"/^" dt (3) 

of a general conformal field theory in two-dimensions. This was achieved by finding a closed form expression for the 
generating function of the moments (p") of p, from which the probability distribution was obtained by inverting a 
Laplace transform. The resulting distribution is conveniently expressed in terms of the dimensionless variable x = pr^ 
and is a shifted Gamma distribution: 

P{x) = ^{x + xo) ^ ^ 1'^ exp(-/?(x + xo)), (4) 



with parameters 



r(a) 



^0=^, a^^, f3^n. (5) 



Here x = —xq is the infimum of the support of the probability distribution, which we will often call the lower bound 
of the distribution, and c > is the central charge, which is equal to unity for the massless scalar field. Using the 
binomial theorem and standard integrals, the n'th moment 

On = / P{x) dx , (6) 



of P is easily found to be 

= E [k) nk + a) = i-xor 2Fo(«, -n- (/Ja^o)"^), (7) 

where 2F0 is a generalized hypergeometric function. 

The lower bound, —xq, for the probability distribution for energy density fluctuations in the vacuum for c = 1 is 
exactly Flanagan's optimum lower bound, Eq. ([2]), on the Gaussian sampled expectation value and, for all c > 0, 
coincides with the result of Ref. [26 . As was argued in Ref. [T5], this is a general feature, giving a deep connection 
between quantum inequality bounds and stress tensor probability distributions. The quantum inequality bound is 
the lowest eigenvalue of the sampled operator, and is hence the lowest possible expectation value and the smallest 
result which can be found in a measurement. That the probability distribution for vacuum fluctuations actually 
extends down to this value is more subtle and depends upon special properties of the vacuum state. In essence, the 
Reeh-Schlieder theorem implies a nonzero overlap between the vacuum and the generalized eigenstate of the sampled 
operator with the lowest eigenvalue. 

There is no upper bound on the support of P{x), as arbitrarily large values of the energy density can arise in 
vacuum fluctuations. Nonetheless, for the massless scalar field, negative values are much more likely; 84% of the time, 
a measurement of the Gaussian averaged energy density will produce a negative value. However, the positive values 
found the remaining 16% of the time will typically be much larger, and the average [first moment of P{x)] will be 
zero. 

The asymptotic positive tail of P{x) has recently been used by Carlip et al [501 to draw conclusions about the 
small scale structure of spacetime in a two-dimensional model. These authors argue that large positive energy density 
fluctuations tend to focus light rays on small scales, and cause spacetime to break into many causally disconnected 
domains at scales somewhat above the Planck length. 

In Ref. |15j . we also reported on calculations of the moments of '.if^: averaged with a Lorentzian, where Lp is & 
massless scalar field in four-dimensional spacetime. It appears that the probability distribution is also a shifted 
gamma function in this case. Define a dimensionless variable x by 

1*00 

X={47TTf f{t)^'dt, (8) 
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where 



There is good evidence that the probabihty distribution is to be Eq. Q with the parameters 

a=4 ^"-^ (10) 

These parameters were determined empirically by fitting to the first three calculated moments. However, the resulting 
distribution matches the first sixty-five moments exactly (agreement had been checked up to the twentieth moment 
at the time of writing of Ref. [E]), so there can little doubt that it is correct. The details of this calculation are given 
in Sect. |III| and Appendix [5] 

Furthermore, the probability distribution for both the two-dimensional stress tensor and the four-dimensional :(/?^: 
is uniquely determined by its moments, as a consequence of the Hamburger moment theorem [H]. This states that 
if a„ is the n-th moment of a probability distribution P{x), then there is no other probability distribution with the 
same moments provided there exist constants C and D such that 

\an\<CD''n\ (11) 

for all n. This condition is a sufficient, although not necessary, condition for uniqueness, and is fulfilled by the 
moments of the shifted Gamma distribution. The Hamburger moment theorem is also an existence result: given a 
real sequence {a„}, n = 0, 1, 2, • • • with oq = 1, such that the N x TV-matrix Hmn = am+n (0 < m, n < — 1) is 
strictly positive definite for every N = 1,2, . . ., then there exists at least one associated probability distribution for 
which the a„ are the moments. 

III. MOMENTS AND MOMENT GENERATING FUNCTIONS 

A. Explicit Calculation of Moments 

In this section, we describe how the moments of a quadratic quantum operator may be calculated explicitly. Let 
be a free quantum field or a derivative of a free field, and let T be the smeared normal ordered square of </>: 



T = y ■4':{x)f{x)dx, (12) 

where / is a sampling function. In our detailed calculations, the smearing will be in time only, and / = f{t) will be 
the Lorentzian function of Eq. ([9|, but our preliminary discussion can be more general. The 7i-th moment /i„ of T is 
formed by smearing the vacuum expectation value 

Gn[xi, . . . ,x„) = {-.(f^-Xxi) ■ ■ ■ ■.(j>'^:{xn)) (13) 

over n copies of /. By Wick's theorem, this quantity is equal to the sum of all contractions of the form 

I 1 I 1 

(l3{xi)(j){xi)(j}{X2)<l){X2)(l)ix3)(l}{x3) ■ ■ ■ (j){Xn)(p{Xn) ■ (14) 



The contractions are subject to the rules that no 4>{xi) is contracted with the other copy of itself and all fields are 
contracted, with each contraction 



I 1 

(t){xi) ■ ■ ■ 4>{xj) (15) 

contributing a factor {<f>{xi)(j){xj)) . 

It is convenient to represent the contractions by graphs with n vertices labelled Xi, . . . ^Xn placed in order from left 
to right so that (1) every vertex is met by exactly two lines; (2) every line is directed, pointing to the right; (3) no 
vertex is connected to itself by a line. For each graph every line from Xi to Xj contributes the factor {cf){xi)4>{xj)) and 
we supply a combinatorial factor that gives the number of contractions represented by a given graph; we then sum 
over all distinct graphs of the above type to obtain G„(xi, . . . , a;„). For example, the graph in Fig. [T^ describes the 
two contractions which contribute to the second moment 



/i2 = 2 / dxidx2f{xi)J{x2){4>{xi)4){x2)) , (16) 
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FIG. 1: The graphs for n = 2 (a), and n — 3 (b) are illustrated. 



SO the combinatorial factor for n = 2 is 2, while Fig. nb corresponds to the eight contractions pairing a 4'{xi) with a 
4'{x2), a (pixi) with a (j){x3) and a (j>{x2) with a 0(a;3)7e.g., 



^(a:;i)(/>(a;i)^(a;2)</>(a;2)(/'(x3)0(x3) (17) 



and 

I 1 I 1 

4>{xi)(l3{xi)(j){x2)(l)iX2)(l)iX3)(j){X3) . (18) 



Note that the moments /z„ have dimensions of inverse powers of length, which depend upon the specific choice of (/). 
It is convenient to rescale the /i„ and define dimensionless moments a„. Our explicit calculations of moments assume 
the Lorcntzian sampling function of width r given in Eq. ^. In the case that 4> — ip^ the massless scalar field in four 
dimensions, we take 

an = (4^t)2« ^„ . (19) 

For the case that 4> — (f, we take 

a„ = (47rr2)2"^„. (20) 

We also take the latter form for the cases of the squared electric field, and scalar and electromagnetic field energy 
densities. 



B. Moment Generating Functions 



For n > 4, the Wick expansion involves both connected and disconnected graphs. However, we need not consider 
the disconnected graphs explicitly, as the moment generating function M is the exponential of W, the generating 
function for the connected graphs. The full moment generating function is defined by 



= E ^ ' (21) 
SO the n-th moment has the expression 

The connected moment generating function, W, has an analogous definition, but in terms of the dimensionless 
connected moments C„ only. These are the moments which arise from counting only connected graphs. For n — 2, 
there is a single connected graph, with combinatorial factor 1 as already described. For n > 2, there are -^{n — 1)! 
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distinct connected graphs, each with a combinatorial factor 2" [?D]. Of course, the enumeration of these graphs 
becomes rapidly unmanageable, and one must exploit further degeneracies among the graphs to reduce the counting. 
For sampling using the Lorentzian function, it is possible to reduce the number of terms to the number of distinct 
partitions of n into an even number of terms. This grows much more slowly than ^{n — 1)!: for example, for n — 30, 
we require 2811 terms instead of 29!/2 « 4.4 x 10^°. Further details can be found in Appendix [A| 

Our procedure will be to explicitly compute a finite number N of connected moments, which allows W to be 
approximated as an iV-th degree polynomial in A. We then use 

M = e^ (23) 

to find M, which may also be approximated as an A''-th degree polynomial. Finally, the first N moments may 
be read off from the coefficients of this polynomial. We emphasize that this procedure makes sense whether or not 



the series (21 1 converges; expressions such as (23) are simply convenient expressions for the combinatorial relation 
between different moments and may be understood as formal power series. 

Consider the case of in four dimensions, where iy9 is a massless scalar field and the average is in the time direction 
only. The two-point function which appears in the integrals for the moments is now 

1 1 f°° 

mW)) = - - 4^ I ^-e-(-* . (24) 

The corresponding dimensionless moments were calculated using MAPLE for iV < 65, and the resulting moments up 
to = 23 are listed in the first column of Table [l] Our computations were exact and give the a„ as rational numbers. 
However, for ease of display, the results have been rounded to five significant figures. The full set of exact moments is 
available as Supplementary Material '31]. As stated earlier, these moments may be used to infer that the probability 
distribution of the quantity in ([s]) is a shifted gamma given by Eqs. ([4]) and (10 1. Only the first three moments are 



needed for this fit, but the result reproduces the first 65 moments exactly, a spectacular agreement. 

Next we turn to the case where ^ = and calculate several of the moments of the Lorentz-smearing of ^P' . In this 
case we use 

3 1 ("^ 

mm) = kmm) = = rf^a^e-^*-* -) . (25) 

As before, the moments were computed exactly as rational numbers using MAPLE for < 65 [31], and the resulting 
moments up to A^ = 23 are listed in the second column of Table |lj 

Once we have a finite set of moments for t/;^, we can calculate the corresponding moments for several other operators 
of physical interest: we give the examples of the energy densities for the massless scalar and electromagnetic fields, 
and the squares of the electric and magnetic field strengths as particular examples. These all take the form 

/oo 
d</(i)^a/:0?:(t,O), (26) 
-oo J 

where aj are constants and the 0/ are (components of) free fields [in the sense that the Wick expansion is valid] with 
two point functions obeying 

c/(5/j(99(i,x)99(i',x'))x=x'=o (27) 

in the vacuum state, where (p is the massless free scalar field as before and the c/ are constants. Defining the 
dimensionless moments for A in the same way as for (^^, one easily sees that the contribution of any connected 
diagram becomes a sum over / of the contributions from each species </)/, with no cross terms mixing different species 
in any given term. Thus 

C„(A)=^(a/C/)"C„(^2)^ (28) 
/ 

from which we may infer 

WiA, A) = £ ^l^A^ = J2 W{^\a,c,\) (29) 



n=0 

and 



M{A, A) = e'^^^^^) = W M{^^, aiciX) . (30) 
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TABLE I: Lorentzian smearings of the Wick square of the free massless field tp^ , the Wick square of its time derivative <^^, 
the square of the electric field strength E"^ , and the energy densities of the scalar and electromagnetic fields ps and Pem 
respectively. 



n 








ps 


pEM 





1 


1 


1 


1 


1 


1 

















2 


2 


9/2 


6 


3/2 


3 


3 


48 


1890 


1680 


525/2 


420 


4 


1740 


2.5516 X 10'^ 


1.5121 X 10*^ 


1.6538 X 10^ 


1.8903 X lO'' 


5 


83904 


8.5527 X 10^ 


3.3789 X 10^ 


2.7057 X 10* 


2.1119 X 10* 


6 


5.0516 X 10'^ 


6.0498 X 10" 


1.5934 X 10" 


9.4918 X 10" 


4.9794 X 10" 


7 


3.6472 X 10* 


7.9890 X 10^^ 


1.4027 X 10" 


6.2499 X lO^-'^ 


2.1918 X 10^^ 


8 


3.0708 X 10^° 


1.7862 X 10^^ 


2.0908 X 10^^ 


6.9804 X 10^^ 


1.6334 X lO^'-* 


9 


2.9538 X 10^^ 


6.2613 X 10^® 


4.8861 X 10^^ 


1.2231 X 10^* 


1.9086 X 10" 


10 


3.1956 X 10^'' 


3.2427 X 10^^ 


1.6870 X 10^'' 


3.1669 X 10^* 


3.2949 X 10^^ 


11 


3.8406 X 10^** 


2.3696 X lO^'^ 


8.2184 X 10^* 


1.1570 X 10^^ 


8.0257 X 10^^ 


12 


5.0767 X 10^* 


2.3561 X 10*^ 


5.4477 X 10^^ 


5.7522 X 10" 


2.6600 X 10^*^ 


13 


7.3196 X 10^° 


3.0960 X lO'*® 


4.7723 X 10" 


3.7793 X 10''^ 


1.1651 X 10"*^ 


14 


1.1432 X 10^3 


5.2487 X 10^^ 


5.3938 X lO''^ 


3.2036 X lO"*^ 


6.5843 X 10*^ 


15 


1.9226 X 10^^ 


1.1252 X 10" 


7.7085 X 10^* 


3.4338 X 10^^ 


4.7049 X 10^° 


16 


3.4641 X lO^^^ 


2.9981 X 10*^^ 


1.3693 X lO''" 


4.5748 X 10" 


4.1789 X lO'^^ 


17 


6.6572 X 10^9 


9.7841 X lO'^^ 


2.9791 X 10*^^ 


7.4647 X 10^^ 


4.5458 X 10®' 


18 


1.3592 X IQp^ 


3.8605 X 10" 


7.8364 X 10™ 


1.4726 X 10'^* 


5.9787 X lO''^ 


19 


2.9384 X 10=*'' 


1.8209 X 10^^ 


2.4642 X 10™ 


3.4730 X 10" 


9.4000 X 10™ 


20 


6.7046 X 10^'' 


1.0164 X lO'*'^ 


9.1702 X 10**^ 


9.6935 X 10" 


1.7491 X 10™ 


21 


1.6103 X 10^^ 


6.6549 X lO''" 


4.0026 X 10*^ 


3.1733 X 10** 


3.8172 X 10*^ 


22 


4.0607 X 10''^ 


5.0695 X lO'*® 


2.0327 X 10^^ 


1.2087 X lO''" 


9.6927 X 10*'^ 


23 


1.0727 X 10'*'' 


4.4604 X 10"^ 


1.1923 X lO^'-* 


5.3172 X 10^^ 


2.8427 X 10^^ 



These results hold for arbitrary smearing on the time axis. 

This procedure may be apphed to the energy density operator for the massless scalar field 



because 



(^(t)aV(i'))x=x'=o = 



which is seen by direct computation of the left-hand side and comparison with Eq. ( 25 1 . Thus we find 



the factor of 3 appearing in one of the numerators corresponds to the spatial dimension. Thus 

W{ps, X)^W (^ip^ ^X^ + 3W (^ip^ 

and 



(31) 

(32) 
(33) 

(34) 
(35) 



M{ps,X)^M{ip',-X) 



M [ -X 



(36) 
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Again, these results should be understood as a relation between formal power series. Concretely, given the first N 
moments of <^^, we can approximate M{(fP',X) as a polynomial, and then use the above relation to find the first N 
moments of ps- The results are tabulated in the fourth column of Table |l] 

Similarly, the components of the square of the electric and magnetic field strength Ei and Bi obey 



{E,{t)E,{t'))^-. 
Following the same line of reasoning as before, we find 



- {B,{t) Bj(i')>x=x'=o - 3 5,, ((^(t)<^(t'))x=x'=o • 



and 



WiE"^, A) = W{B'^,\) = iW -A 



M{E^,\) = M{B'^,X) = 



M ip^, -A 



(37) 

(38) 
(39) 

(40) 



This result leads to the moments of the squared electric field, tabulated in the third column in Table |lj The results 
for the square of the magnetic field are identical. 
Finally, because we also have 



(S,;(t,x)S,(t',x'))x=x'=0 = 0, 

the energy density of the electromagnetic field 

1 



has connected moments 



and hence 



CniPEAl) 



Pem 



{E^ + B^ 



CniE^ 



1 



WipEM, X) = 2W[ E\ ^A ) = ( \\ 



and 



M{pEM.X) 



m[e\-x 



M(^2 -A 



(41) 
(42) 

(43) 
(44) 
(45) 



leading to the remaining entries in Table |T] 

An important observation is that these moments (apart from those of the Wick square) grow too rapidly to satisfy 
the Hamburger moment criterion, Eq. . This may be confirmed by noting that in all cases In a„ grows faster with 
increasing n than n Inn + CiTZ + Cq for any constants Cq and Ci. In fact, the growth for ip^ is shown in Appendix [b| to 
be of the form 



CD" (3n-4)!, 



(46) 



where the constant D is proved to lie in the range 3.221667 < D < 3.616898 (our numerical evidence suggests 
D ~ 3.3586). For probability distributions known to be confined to a half-line, which is the case here, there is 
a sufficient condition for uniqueness which is weaker than the Hamburger moment criterion. This is the Stieltjes 
criterion [19|, which is 

a„ < CL>"(2n)!. (47) 
Unfortunately, this criterion is also not fulfilled here. This means that we cannot be guaranteed of finding a unique 



probability distribution P{x) from these moments. This issue will be discussed further in Sec. VHIA 

Note that in four dimensions, the operators {ip'^, E^, ps, and pem) all have dimensions of length^"^. Their 
probability distributions P{x) will be taken to be functions of the dimensionless variable [See Eq. (20 1.] 

X = (477X^)2 A, 

where A is the Lorentzian time average of (tp'^, E^, ps, Pem)- 



(48) 
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C. Lower Bounds 



In general, we may use relations between different moment generating functions to find relations between the 
corresponding probability distributions, and especially between the lower bounds of these distributions. (Strictly, 
these are the infima of the support of the distributions.) Let p{x) and q{x) be two probability distributions, with 
moment generating functions M{p, A) and M{q, A), respectively. These generating functions can be expressed in terms 
of the bilateral Laplace transforms of their probability distributions: 



J 



and 



M{q, A) 



q{x) e dx 



(49) 



(50) 



These integrals are guaranteed to converge at the lower limits, due to the lower bounds on the support of our 
probability distributions. To assure convergence at the upper limit, we may assume Re A < 0. However, many of our 
arguments below do not require convergence of the integrals, which may be regarded as formal power series in A on 
replacing the exponential by its Taylor series. Now let p'¥q[x) be a probability distribution defined as the convolution 
of p and q: 



p * q{x) 



dx' p{x — x') q{x') 



(51) 



As is well-known in probability theory, this is the distribution for the random variable obtained as the sum of 
independent random variables with distributions p and q, and its moment generating function is 



M(p * g. A) — 



dx 
dx' 



dx'p{x~x')q{x')e'''= 



dxp{x-x')e^^'"'''^ 



q{x')e^^ 



dup{u) 



dx q{x') e 



M(p, A)M(g,A) 



(52) 



where u = x — x' . Thus the moment generating function of a convolution is the product of the individual generating 
functions; again, this holds in the sense of formal power series, irrespective of convergence issues. 

We can also give the relation of the lower bounds. As is also well-known in probability theory, the support of a 
convolution p * q of two distributions consists of all values expressible as the sum of a value in the support of p and 
a value in the support of q. In particular, the greatest lower bound on the support is the sum of the lower bounds of 
the individual distributions. Explicitly, if bp and bq be the lower bounds of p and q, then 



p{x) = if x < bp ; q{x) ~ if x < bq . 
The integrand of Eq. (51 1 vanishes if either x' < bq, oi x — x' < bp. This implies that 

p * q{x) — if X < b.„ 



(53) 



Up V bq . 



(54) 



and in fact this is the greatest lower bound. Thus the lower bound of p * g is the sum of the bounds of p and of q. 
Next consider the effect of a rescaling of A, and let Pa{x) — \a\p(ax), where a 7^ 0. Then 



M{p^,\) = \a\ j p{ax) e^^ ^x ^ J p{x') e^^/")^' dx' = M (p, ^ 



(55) 



Provided a > 0, Pct = if a; < bp/ a, so the effect of rescaling A in M is a rescaling of the lower bound by the same 
factor. If a < 0, the lower bound on the support of Pa is — times the upper bound on the support of p, if this 
exists; if there is no upper bound on the support of p, then evidently pa has no lower bound in this case. 

Now we may combine these results to relate the lower bounds of various probability distributions to that for ip^ . 
Applied to a general operator of the form ( 26 ) , they suggest that the probability distribution for A is a convolution 



of several copies of the probability distribution for ip , with various scalings. For example, Eq. (36) suggests that the 
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probability distribution for the energy density ps, smeared along the time axis, is equal to the convolution of four 
copies of the probability distribution for (^^, with various scalings. In particular, recalling that xq{A) denotes the 
greatest lower bound on the support of the distribution for A smeared in time, this suggests that 

xo{A)=[Y,aici\xQ{^^) (56) 



Hence Eq. (36l suggests that xo{ps) =~ (1/2 + 3 x l/6)xo{(p'^) = xoiip'^). Similarly, Eq. (40) suggests that xo{E'^) = 



3 X (2/3)a;o(^) = 2xo{(p'^), and Eq. (45) suggests that xo{pem) = 2 x {1/2)xo{E^) — x^(E^). In summary, 

xoipEAi) = xo{E^) = 2 xqIps) = 2 xo{(p'^) . (57) 

Likewise, if we consider a combination such as the pressure Tn — \{'f^ + {difY ~ (92^)'^ — (^av)^), we obtain the 
expected result that the probability distribution is unbounded both from above and below. The above derivations 
should be take as suggestive, rather than rigorous proofs, because of concerns about the uniqueness of the underlying 
probability distributions. However, it would be possible to prove them by writing the smeared operator for ps, for 
example, as a sum of mutually commuting self-adjoint operators, each of which was essentially a multiple of the 
smeared operator (under a suitable unitary transformation). This could be done by writing the field in a basis 
of spherical harmonics, in this framework, the three powers of M{(p'^,X/6) arise from the £ = 1 angular momentum 
sector, while the single power of A'I{tp'^, A/2) arises from the £ = sector. Indeed, one of the first quantum inequality 



bounds on the expectation value of ps used precisely this decomposition pS] . More generally, Eq. ( 27 ) could be used 
in conjunction with Wick's theorem to show that timelike smearings of and :0j: commute for / J, at least in 
matrix elements between states obtained from the vacuum by applying polynomials of smeared fields, and might be 
used to put the other relationships above on a firmer footing; we will not pursue this here. 



IV. LOWER BOUND ESTIMATES 



Here we will discuss a technique, a Stieltjes moment test, by which knowledge of a finite number of moments may 
be used to obtain an approximate estimate of the lower bound. If P{x) is a probability distribution with a lower 
bound at a; = —xq, then its moments are 

oo 

a;" P{x) dx . (58) 

-Xo 

Let 

/•OO 

/(y)= / {x + v)\q{x)\^P{x)dx, (59) 



where q{x) is a polynomial and y > xq. We see that I{y) > because the integrand in Eq. (59) is non-negative. If 

N-l 



q{x)=J2PnX\ (60) 



n=0 

then 

N-l 

I{y)^ J2 M^^iN,y)f3*^l3n>0, (61) 

m,n— 

where M{N,y) is a real symmetric N x N matrix with elements 

Mmn{N, y) = am+n+l + ya-m+n {0 < TTl, n < N - I) . (62) 

Let /3„ be the components of an eigenvector with eigenvalue A, then '^^^q Mmn{N , y) Pn = A /3m, and I{y) = 
A J2m=o l/3mP- It follows that M{N, y) has no negative eigenvalues, that is, it is a positive semidefinite matrix, which 
we denote by M{N,y) > 0. This holds for all N and all y > xa. However, as y decreases below xq, the lowest 
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TABLE II: Table of the lower bounds, j/jv, for both tp^ and ip^ 



N 






2 




n ni 071401 940 


3 




ni 41 4954099 


4 


1 94-78'^Q8SfiO 


ni 584995314 


5 




ni 69ni 99565 


6 


0.13872875370 


0.01762865715 


7 


0.14271593142 


0.01816742316 


8 


0.14570717836 


0.01858660399 


9 


0.14803421582 


0.01892432539 


10 


0.14989616852 


0.01920370321 


11 


0.15141979779 


0.01943965011 


12 


0.15268963564 


0.01964226267 



TV 


yN{^ ) 


yN\^ ) 


13 


n 1 5376491 Sn5 


n 01 981 864633 


14 


n 1 5468536476 


n ni 997396948 


15 


n 1 5548374879 


n n9ni 1 9n6075 


16 


n 1 561 8937796 


n n9n93587746 


17 


0.15679884907 


0.02034769569 


18 


0.15734684979 


0.02044932047 


19 


0.15783718730 


0.02054219985 


20 


0.15827850807 


0.02062751059 


21 


0.15867781217 


0.02070622001 


22 


0.15904082736 


0.02077913144 


23 


0.15937228553 


0.02084691828 



N 


11 TIT { J 


7(Ar ( 

yN\^ ) 


24 


n 1 596761 301 8 


09091 01 4970 


25 


n 1 59955674nn 


09096931 050 

KJ .\J^KJlJ\JC/0 ±\JO\J 


26 


n 1 6n91 379090 


091 09481 644 

\J .\J ^ ±\JZj^iJ ± \J ± A. 


27 


0.16045265677 


0.02107702642 


28 


16067453067 


0.02112625203 


29 


0.16088110659 


0.02117276528 


30 


0.16107391397 


0.02121680481 


31 


0.16125428495 


0.02125858099 


32 


0.16142338519 


0.02129828002 









eigenvalue is eventually zero and then negative eigenvalues can occur. Define as the minimum value of y at which 
M{N,y) > 0; in practice, it is easiest to compute y^ as the largest root of the A^'th degree polynomial equation 

detM(iV, y) = . (63) 

Because M(N, y) is a leading principal minor of M (A^+l, y), M{N+1, y) > implies that M{N, y) > 0. Consequently, 
Vn+i > UN and the sequence in N converges to a limit with 

?/oo = hm yN <xo- (64) 

AT— >-oo 

Given a set of moments a„ of an unknown probability distribution, we may form the matrices M{N,y) as above 
and determine the values of j/at. The above argument shows that if j/^v — >■ oo then the a„ cannot be the moments of 
a probability distribution whose support is bounded from below. On the other hand, suppose that a finite limit j/oo 
exists. Then for any probability distribution P with the same moments and support bounded below by —xq, we have 
Vac ^ io- In particular, there is no probability distribution accounting for the given moments with support contained 
in (-?/oo,oo). 



Let us first apply this method to the case of the ip'^ distribution, given by Eqs. Q and (10 1, for which the exact 
lower bound is known. The results of the calculation of the y^ through = 32 are given in Table [Tl| (computations 
were performed in MAPLE to 40 digit accuracy; the reported rounded figures are stable under increase of the 
number of digits). We can improve the estimate of the lower bound by extrapolation. A trial function of the form 
yN = a + b/N + c/N"^ and a least-squares fit using MAPLE |1T] to determine values of a, b and c leads to 

, 2^ 0.167821368174 0.001164170336 , , 

yNiv ) ~ 0.166666666057 + —3 (65) 

The above fit was obtained using the data points for 21 < TV < 33, with residuals of order lO"^'^ over these values, 
and no more than 1.1 X 10"*^ for 2 < iV < 20. Using the fit displayed above, our lower bound estimate now becomes 



Voo = 0.166666666057, in extremely good agreement with the exact bound, xq = 1/6, obtained from Eq. (10 1. This 
suggests the conjecture that —yoo might also coincide with the lower bound of the probability distribution in other 
cases as well, but note the caveat at the end of this section. 

A different numerical approach is to use an accelerated convergence trick: given any sequence y — (yN), define a 
new sequence L^'^^y with terms 

{L^''^y)N = iVN+i - Vn) + VN] (66) 

for finite sequences, L'^^'^y is one term shorter than y. This is a linear map on sequences, preserving constants and 
acting onyiq — l/N^ by 



(67) 
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TABLE III: Table of the accelerated lower bounds for both ip^ and ip^ . 



N 




r(3/2)r(l)r(l/2) / -2^, 


21 


n 1 fififififi"! 39^14- 


n9'^(S1 4791 9"^ 


22 


n 1 fififififi^^fii 1 

W.-LL'UU'JUUUU-L-L 


fl 09361 451 051 


23 


n 1 fififififi'ifiQQS 


fl fl93fi1 439088 


24 


n 1 fififififi'iSI 

U. -LULUJUvJUOJ-tJO 


fl fl93fi1414Q7S 


25 


n 1 fififififiSQl 35 


fl fl93fi1 399500 


26 


0.16666659972 


0.02361385460 


27 


0.16666660689 


0.02361372693 


28 


0.16666661307 


0.02361361053 


29 


0.16666661843 


0.02361350416 


30 


0.16666662310 


0.02361340672 


31 


0.16666662718 





for any p, fc > 0. Thus if un = a+bN^''+cN~^+- ■ ■ , with £ > k, the sequence L'^'^^y converges to a as 0{N^ min{£,fc+i}-j^ 
rather than 0{N~''). This trick may be repeated: in the situation above, L^^^L*^^)y(</J^)jv would be expected to 
converge with 0{N^^) speed to the limit. The results give values differing from 1/6 by less than 10~^ for all 
11 < < 31. Part of the 'accelerated' sequence is given in Tabic [Till 

We may now apply the same procedure to the case of (p^, where the exact bound is not known. The yN^'fi'^) are 
also given in Table pi} and clearly converge more slowly than those of the y^if'^)- Indeed, successive differences 
VN+ii^"^) — l/Ni'-p'^) appear to decay as 0{N~^/^). A least squares fit to the trial function yN{v^) = a + h/N^/"^ + 
c/N + d/7V3/2 

gives 

, 2, 0.012425890959 0.002768353926 0.006533917931 , , 

y^i^') « 0.0236174942666 (68) 

using 21 < iV < 33, with residuals less than 1.2 x 10^^" on these values, and no more than 10^^ on 6 < < 20. 
Applying the acceleration technique, L'-^^^^ L'-^^ L'^^^^^ y{ip'^) gives a sequence differing from 0.02361 by no more than 
8.1 X 10"'^ on 11 < TV < 30. Taking this together with the least squares fit gives reasonable confidence in an estimate 
y^{ip^) = 0.02361 ± 1 X 10^^ 

In contrast, the non-optimal bound for (p'^ and ps, given by the method of Fewster and Eveson [27j, is xo{FE) = 
27/128 « 0.21, which is an order of magnitude larger. [This bound is given by minus the right hand side of Eq. (5.6) 
in Ref. [57] multiplied by (47rr^)^.] If, in fact, —j/oo coincides with the lower bound of the probability distribution, we 



can now use the results in Eq. (57 1 to write our estimates of the probability distribution lower bounds as 



- xo{pem) = -xo{E^) « -0.0472 - xo{ps) = -xo{(p^) « -0.0236 . (69) 

These are also estimates of the optimal quantum inequality bounds for each field. 

There is an important caveat to this reasoning, however. If the moments do not correspond to a unique probability 
distribution (i.e., if it they are indeterminate in the Hamburger sense) then there will exist probability distributions, 
called von Neumann solutions in Ref. |19j . with the given moments that are pure point measures, in contrast to the 
continuum probability distribution that would be expected for the quantum field theory operators we study (and 
which we find for the c/)^ case). As the moments arise from a probability distribution supported in a half- line, there 
is a distinguished von Neumann solution, called the Friedrichs solution in Ref. |19| . that is supported in a half- line 
[—xf,oo) and has the property that no other solution to the moment problem can also be supported in [—xf,oo). 
(See Appendix CI of "19] for a brief summary.) Hence if operator A has Hamburger-indeterminate moments, we would 



have yoo{A) = xf{A) < xo(A). Nonetheless, the results in Eq. (69) would still be true if the approximation signs are 
replaced by <. 



It is of interest to note that the magnitudes of the dimensionless lower bounds, given in Eq. (69) are small compared 
to unity. Given that the probability distribution must have a unit zcroth moment and a vanishing first moment, this 
implies that P{x) ^ 1 in at least part of the interval —xq < x < 0. Thus the spike at the lower bound found 
in the two-dimensional case may be a generic feature. The small magnitudes of xoifs) a-nd xo(pem) imply strong 
constraints on the magnitude of negative energy which can arise either as an expectation value in an arbitrary state, 
or as a fluctuation in the vacuum. They also imply that an individual measurement of the sampled energy density in 
the vacuum state is very likely to yield a negative value. 
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FITS FOR THE APPROXIMATE FORM OF THE PROBABILITY DISTRIBUTION 



In this section, we explore the extent to which knowledge of a finite set of moments may be used to obtain information 
about P{x) beyond the lower bounds found in Sec. IV 



A. A procedure to find the parameters of the tail of P{x) 

We begin with the large x limit. Let us adopt the ansatz that: 

P(a;) -coa;^e-°"% (70) 
for large x. We assume that we can use this form of the tail to compute the large n moments, and find 

P{x) dx w Co 

= ^a-("+''+i)/^ [(n + & + l)/c- 1]!, (71) 

for n ^ 1. We expect the dominant contribution to come from a; ^ 1, so we set the lower limit in the second integral 
to zero for convenience. 



Next we compare Eq. (71) with the Eq. (461 for the large n form of the moments. This comparison reveals that 
we should have 



-2 a^D"^/^ co^CD/S. 



With these values for c and b, the ratio of successive moments from Eq. (71 1 becomes 

a„+i 3(n — l)(3n — 2)(3n — 1) 



(72) 



(73) 



Now we may use the computed values of two successive moments, such as n = 64 and ti = 65, to find the value of 
a, and then the value of cq from Eq. (71). The results for the different operators are listed in Table IV It should 



be noted that knowledge for further moments beyond n = 65 could change the values in this table. A rough error 
analysis suggests that these values are correct to about five significant figures. 



TABLE IV: Values of the Parameters for the Tails, in the form of Eq. (701 



Operator 


Co 


a 


b 


c 




0.47769605 


0.6677494904 


-2 


1/3 




0.95539211 


0.7643823521 


-2 


1/3 


PS 


0.23884802 


0.8413116390 


-2 


1/3 


Pem 


0.95539211 


0.9630614156 


-2 


1/3 



The values of the constants a and cq for the various cases can be related to one another by means of the relations 
between the connected moments, Eqs. (34 1 and (43), derived in Sect. IIIB First, we need the fact that the connected 



moments and the full moments rapidly approach one another for large n, specifically 

C„-a„(l + 0(7i-4)), 71 >1. 



(74) 



This relation may be demonstrated analytically, or inferred numerically from the computed moments. This means 
that Eqs. (34) and (43) hold for the full moments, a„ when n is large. The former relation may be simplified to 



o-n{ps) ~ 2 "a„(95 ). The asymptotic form, Eq. (46), for the moments of if also holds for the other operators, but 
with different choices of the constants C and D: 



C((^2) ^ c{ps) = ^CiE') = ^C{pem) and D{^') = 2D{ps) = ^DiE') = 3D{pem) ■ 



(75) 
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For example, anipEAi) = 6 (1/3)" a„((/3^), from Eq. (43), implies the above relations between C{ip^) and C{pem) and 
between D{ip'^) and D{pem)- These relations and Eq. (72) imply that 



co{ip^) = 2co(ps) = ^co(S^) = ]^co{pem) and a((^^) = 2 3 a{ps) 



aiE'')^3--^aipEM). (76) 



These relations are borne out by the values in Table [TV] 



B. Estimating when our tail fit is a good approximation 



Since the Hamburger and Stieltjes moment conditions are not fulfilled for our moments, we do not know whether 
our probability distributions are unique. However, if we assume that they are, then we can estimate the range in x 
where we expect our fitted tails to give a good estimate of the actual distributions. Our general form for the tails of 
the probability distributions is approximately 

Pfit(x)^coa;-2e— (77) 

As an example, for pem, this gives a good fit (< 10%) for n = 4, 5, 6, 7, 8 and a better fit (< 1%) for 9 < n < 64. (We 
used n = 65 to set cq, so it should not count.) Let 

/„ (x) = a;" Pfit (a;) = Co e— (78) 



so 



fn{x)dx 



is our predicted moment from the above form. The maximum of the function fn{x) will be where fni^) 
corresponding to 

'3(n- 2) 



(79) 
= 0, 

(80) 



If Pmix) gives a good approximation for An, then it should give a good approximation to the exact P{x) for x ~ 

^(•^max) ■ 

For the electromagnetic energy density a « 1, so for n — 4, a;niax ~ 216, and for n = 65, Xmax ~ 6751269. Thus 
if Pfit gives reasonable fits to the moments for 4 < n < 65, then it should be a fair approximation to the exact 
distribution in the range, roughly, 10^ < a; < 10^, assuming uniqueness of the distribution. 



C. Approximate fits for P{x) including the inner part 

One can attempt to model the entire probability distribution, including the inner part, by experimenting with 
functions of the form: 

P{x) ^ ci {xq + x) " exp[— /3(xo + x)''] + Co (ao + (xo + 2;)^) exp[— a(xo + a;)^^^] . (81) 

A reason for using this form is that one need not bother with trying to match inner and outer parts of the function. 
Depending on the choices of the constants, one can possibly get the first term to dominate for small x, and the second 
for large x. We use the values of a, 61 from the tail fits and the values of xq from the quantum inequality bounds 
given earlier in Sec. |IV| 

The most interesting case is the distribution for pem, the electromagnetic energy density. For the values of the 
constants given in Table |vj the fractional errors between the calculated and fitted moments in the O"^ through 22"^^ 
moments are given in Table [ylj Since the exact value of the first moment is 0, we list the fitted value separately as: 
moment = 0.0247001. The errors in the fourth and fifth moments are somewhat large 15%), but the errors 
tend to progressively decrease as we go to large n. So this heuristic model distribution gives a reasonably good fit for 
the innermost part of the distribution and the tail, but does somewhat poorly for the middle part of the distribution. 
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FIG. 2: The graph of P{x) vs x of our fit to the probabihty distribution function for pem, the electromagnetic energy density 
sampled in time with a Lorentzian of width r. Here x = IGtt^t'^ pem ■ The distribution has an integrable singularity at the 
conjectured optimal quantum inequality bound x = —xq — —0.0472. 



The graph of P{x) vs x for this case is given in Fig. [2] It has a spike (an integrable singularity) at the quantum 
inequality lower bound. However, our method may not be sufficiently sensitive to conclude the existence of this 
singularity. It is possible that there are non-singular distributions which fit the first several moments as well as does 
our postulated form. Thus we cannot conclude whether the actual distribution has a spike in it at the lower quantum 
inequality bound, as indicated in the plot. The distributions for ps and (p'^ in two-dimensional spacetime, which are 
known exactly and uniquely, both have a spike at the quantum inequality lower bound, as does the distribution for 
ip'^ in four dimensions [B] . In Tables |V] and |Vl[ we list the fitting constants and fractional errors, respectively, for the 
Pem probability distribution. The values of the constants were obtained by calculating the moments from Eq. (81) 
and using the MATHEMATICA Manipulate command to adjust the values of the constants to get the smallest 
fractional errors between the fitted moments and the actual moments. 



TABLE V: Fitting Constants for the Model Distribution for psM in Eq. (811 



Constant 


Pem 


a 


0.9630614156 


Co 


0.95539211 


Xo 


0.0472 


ceo 


610 


Cl 


0.028 


/3 


19.65 


7 


1.05 


a 


0.9999 



VI. BOUNDS ON THE CUMULATIVE DISTRIBUTION FUNCTION 



As already mentioned, it is possible that the moment problem is indeterminate and that there are many probability 
distributions with these moments. Here, we show that no such distribution can have a tail decreasing much more 
slowly than that studied above. Our tool for this purpose is a simple variant of Chebyshev's inequality: if X is any 
random variable taking values in [— xq, cxo), with moments an, then the probability Prob(X > A) that X exceeds any 
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TABLE VI: Table of Fractional Errors. Here the fractional error is [a„(fit) — a„]/a,i, where the a„ are given in Table and 
the a„(fit) are computed from Eq. (811. For the n = 1 case, the fractional error is not defined, since the first moment is 0. 
Fractional errors in succeeding moments beyond n — 5 are progressively smaller. Although all moments through n = 65 were 
used, we display the fractional errors through n = 21. 



n 


Pem 





0.00450644 


1°* moment 


not applicable 


2 


-0.00661559 


3 


-0.0770297 


4 


-0.152164 


5 


-0.150279 


6 


-0.117773 


7 


-0.0843077 


8 


-0.0590582 


9 


-0.0420107 


10 


-0.0308225 


11 


-0.0233756 


12 


-0.0182526 


13 


-0.0145945 


14 


-0.0118911 


15 


-0.00983456 


16 


-0.0082327 


17 


-0.00696063 


18 


-0.00593416 


19 


-0.00509465 


20 


-0.00440012 


21 


-0.00381978 



given A is bounded by 

Prob(X > A) < "» + P'-"M^<OK (82) 
for all n. To prove this, let d^(x) be the probability measure of X and then compute 

poo pOO pOC pO 

A"Prob(X > A) ^ A" / dfi{x) < / ^"^^(x) < / x"d/.t(x) = a„ - / x"dn{x) < + Prob(X < 0)a;o . (83) 

J\ J\ Jo J -xq 

(The term Prob(X < 0)xq is only needed for odd n, in fact. We have also written d^{x), rather than P{x)dx for the 
probability measure to emphasise that we are not assuming a continuous probabihty density function.) In our case, 
we know that < xq{FE) < 1, so we have 

Prob(X > A) < inf ^!L±I. (84) 

riGN A" 

Now, for moments growing as a„ ^ CI?"(3n — 4)!, the ratio of successive terms in the infimum is 

an+i + 1 ^ (3n-l)(3n-2)(3?i~3) 
X(an + 1) A 

so, for each fixed A, the sequence will decrease until the term where n ^ ^{X/Oy^^ and will increase thereafter. This 
gives an asymptotic bound on the tail probability 

Prob(X> A) <c(^j' r((A/Z?)i/3-3) ^ V2^cf e-(^/^)'''. (86) 
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as A — >■ oo. 

This gives an upper bound on the tail probabihty distribution, which is not much more slowly decaying than that 
for our fitted tail, for which the tail probability would be decaying like C (x)*^"^ g^i^^/D) _ ^j-^g following discussion 
sketches how information on the lower bound can be obtained; this could be developed into a rigorous discussion (and 
probably sharpened) with further work. In fact, we do not seek a strict lower bound on the tail probability, but rather 
aim to show that it must be very often of the order of the fitted tail or higher. 

Let Q[x) = Prob(X > x). Then we have, for any A > A > xo, 

/"OO 

a„ < A"Prob(X < ^) + y x"dfi{x) (87) 
< A"Prob(X < A) + Q(A)A" + n j Q{x)x"-^ dx 

X 

D 



<y' + nl Qix)x''~Ux + V2^CD''-^n I (|- ] " da; (89) 



<X'' + nJ^ Q{x)x''-^ dx + 3\/2^CD"nr(37i - 7/2, {A/DY^^) (90) 

in which we have integrated by parts in the second line and used the fact that Q{X) = 1 — Prob(X < A), as well as 
the upper bound found above; r(A^, z) is the upper incomplete P-function. We can now make n-dependent choices 
of A and A so that the first and third terms are negligible in comparison with a„ for large enough n. For example, 
A = {4:n)^D and A — n^D will do: it is a simple application of Stirling's formula to see that A"/(D"(3n — 4)!) ~ 
const X 7i^/^(e/3)^" — 0; for the upper end we first estimate r(3n — 7/2, 4n) ^ 4(4n)'^"~^/^e~^" using Laplace's 
method (see [32j, section 4.3) [4^ which gives 

_^r(3„ 7/2 4„)^ 1 /3y/'/jiV^0. (91) 



(3n-4)! 72^ V4/ V27e 

With these choices of A and A in force, we set F{x) = xQ{x)e^^^^^^^^ , whereupon we have 

i^(x)x"-2e-(^/-D)'^' dx > CD"(3n - 4)! (92) 



from ([90]). Now let S be the subset of x e [A, A] for which F{x) > \CD{D/xY/^ . We bound F from above by 
y/2^CD{D/xf/^ on S, and by \CD{D/xYl^ on the complement S" of S in [A, A], to give 



i^(a;).x"-2e-(-/^)''' dx < ^CD^^ f ^n-iye^-i./D)-'^ ^ CD^ j ^„_7/3p-(./D)i/^ (93) 

3'= 



Now the first integral on the right-hand side can be bounded from above by the supremum of the integrand multiplied 
by the Lebesgue measure \S\ of S, while the second is bounded by the integral over all [0,cxd). The supremum 
mentioned occurs for x = (3n — 13/2)^£', and we find 

CD"{37i - 4)! < \S\V2TTCD"~\i{3n - i3/2)3"-i3/2e-(3"-i3/2) + icr»"3nP(3n - 4). (94) 

Rearranging and using Stirling's formula, this requires 

(',3i? — 4Vp3"-i3/2 n 27 
' y8^7i(3n-13/2)3»-i3/2 2 ^ ' 

Summarizing, we have shown that in the interval [n^D,An^D], for n sufficiently large, we have 



1 „ /D 



4/3 



Prob(A: >x)>^C i^^j e-^^^l^)"' (96) 

on a set with measure at least "^D-n? . It seems likely that this is a substantial underestimate of the measure of S, as 
some of the estimates used in the last part of the argument are rather weak. 

Thus the broad behavior of the tail of the probability distribution is determined by the moments, even if the exact 
probability distribution is not uniquely determined. In the applications we give below, it is only the broad behavior 
that is required. 
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VII. POSSIBLE APPLICATIONS FOR THE TAIL 



Black Hole Nucleation 



The fact that the energy density probabiUty distribution has a long positive tail implies a finite probability for 
the nucleation of black holes out of the Minkowski vacuum via large, though infrequent positive fluctuations. This 
probability cannot be too large, of course, or it will conflict with observation. Let us sample a spacetime region (a 
cell) over a size ^ « r, where r equals the sampling time. For an energy density p, which is roughly constant in space, 
the associated mass will be M w pi^. This can be a black hole if GM « £, or ip^M w £, in units where h — c — 1 and 
£p is the Planck length, which implies p « l/(£p^£^). Here we chose r so that the sampling time is approximately 
the light travel time across the black hole. 

Note that we should really use the probability distribution for energy density sampled over a spacetime volume, 
with the spatial and temporal dimensions approximately equal. For the purpose of an order of magnitude estimate, 
we assume that the probability distribution for sampling in time alone will yield roughly similar results. 

Let our observation volume be V and our total observation time be T. The number of cells in this spacetime volume 
is iV = VT/£'^. Because black hole nucleation will be a rare event, we assume that different nucleation events will be 
widely separated and uncorrelated. The number of black holes, n, nucleated in this spacetime volume, VT is then 
n w NPn, where P„ is the probability of a black hole nucleation in our sampled spacetime volume i"^. Let us estimate 
that 



Pn 



2x 



P{y)dy. 



(97) 



where 



X = 167r T p = IOtt 



167r^ 



M 



(98) 



and nip is the Planck mass. Here P„ is the probability of nucleating a black hole in the range between x and 2x. 
However, in the limit of large x, Pn will be independent of the exact upper limit in Eq. (97). Let the probability 
distribution have a tail of the form given by Eq. (77). Then 



i'2x /•U2 

Pn~CQ y^^e^^y dy = icaa^ I w~''e""rfw = 3 Cq a^[r(-3, Mi) - r(-3, ua)] ■ 



(99) 



Here u — ay^^"^, Ui = aa;^/'^, U2 = 2^/^ ui, and r(— 3, u) is an incomplete gamma function. This function has the 
asymptotic form 

r(-3,u) ww-^e"" (100) 
for u ^ 1. From this form, we see that the contribution from the lower integration limit dominates, and we have 



Pn 



3 Co 



X 3 e 



(101) 



for large x. 

Thus we have for the mean number of nucleated black holes 



or, using Eq. ( 101 ), 



£4 



ip^M^ 



20/3 



(102) 



(103) 



1/3 

where ao = (IGtt^) a. For the energy density of the EM field, cq ~ 0.955, a « 0.963, so oq « 5.2. Therefore for this 
case we have 



n w 10" 




20/3 



exp[-5.2(Af/mp)^/^] 



(104) 
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To estimate the probability of black hole nucleation, let us first choose V = Icm'^, T = 1 sec, and n = 1, which 
gives VT/ip'^ w (10'^^)'^ 10''"^ w 10^^^. We want the probability of one black hole forming in one cubic centimeter of 



space over an observation time of one second. We can use Eq. (104) to determine the resulting mass of the black 



hole, which turns out to be M w 400 m^. Let us now consider our observation volume and time to be the size and 



age of the universe, which gives VT/£p'^ sa (10^^/10 ^^)^ « 10^''^. Taking n = 1 again, and using Eq. (104 1, yields 



M w 990 TOp. Therefore, if we observe a volume the size of the universe for a time equal to the age of the universe, 
we are likely to see the nucleation of only about one lO^Wp black hole from the vacuum. 

Thus nucleation of black holes of mass ^ 10^ nip is common, but lO'^m^ black holes are very rare. Why do we not 
notice these 400 Wp « 10~^ g black holes? Presumably they must appear for a very short time and be surrounded by 
negative energy which quickly destroys them. 

B. Boltzmann brains 

Recently, the "Boltzmann brain" problem has become the subject of increasing interest in cosmology [331 134) . 
This is the possibility that conscious entities, which may or may not resemble biological brains, might spontaneously 
nucleate and exist for a finite time. Anthropic reasoning requires a count of observers, as the anthropic prediction 
for the value of an observable is the value most likely to be found by a typical observer. If the typical observer is a 
Boltzmann brain in intergalactic space, and not an observer on an earthlike planet, this would greatly alter anthropic 
predictions. As a somewhat more speculative application, we consider what the tails of our probability distributions 
have to say about the probability of nucleating Boltzmann brains in four-dimensional Minkowski spacetime. This 
calculation is similar to the one above for the nucleation of black holes. 

Consider a spatial region of size £, a timescale r, and a mass M, so that the mean energy density is p w M/i^. 
We want to use the tail of the EM energy density probability distribution to estimate the probability of mass M 
appearing in this specific region in a particular interval t. Our sampled energy density is a; = Wtt^t'^ p ~ t'^ M/£^. 
So we have that 

Fix) oc e-^^''' « e^-'^' « exp - r t (105) 

where we have ignored the prefactor and used a ~ 1. The prefactor would contain information about the fraction of 
mass M's that could think. Even if very small, this factor is likely to pale in importance compared to the exponential 
factor derived below. Let M = 1kg « 10^^ cm"-'^, £ = 10 cm, and r — 0.3 sec « 10^" cm. These values give 

^1^^1026, (106) 

so 

P « . (107) 

This is much larger than the exp(— 10^") estimate of Page [3S], who assumes that P oc , where / — Mt = action. 
So our energy density probability distribution makes the Boltzmann brain problem worse. Although the probability 
per unit volume for the nucleation of a Boltzmann brain may seem exceedingly low, the available volume could 
make them more numerous than other observers. Note that in this case, the energy density has been averaged over a 
spacetime region which is much larger in the time direction than in the spatial directions, t ^ £. Hence the probability 
distribution for the energy density averaged in time alone should be a good approximation here. 

VIII. DISCUSSION 
A. Uniqueness Issues 



As was noted in Sec. IIIB[ the moments which we calculate for ■.ip'^: and related operators satisfy neither the 



Hamburger condition, Eq. (Ill, nor the Sticltjes condition, Eq. (47 1 for uniqueness. Thus none of our results for P{x) 
are rigorously guaranteed to be unique. However, there are some observations which are relevant here. First, these 
are sufficient, but not necessary, conditions for uniqueness. There is a necessary and sufficient condition [T^, but this 
condition requires detailed knowledge of all moments and does not seem to be testable in our problem. Second, rapid 
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growth of moments does not automatically mean non-uniqueness. There are examples of sets of moments which grow 
at arbitrary rates, but nonetheless are associated with unique probability distributions. 

On the other hand, if the probability distribution is continuous, with probability density function p{x) on [—Xq, oo), 
and 

\og{p{x))dx ^^^^^ 



then the Stieltjes problem is indeterminate for the moments of p (assuming they all exist, and that Xq < 1 for 
convenience); there is more than one probability distribution supported in [— xq, oo) with the same moments. This is 
a theorem of Krein (modified slightly to our setting) see, e.g., Theorem 5.1 in Ref. |36j . In particular, this would show 

that any distribution whose tail was exactly equal to Pfit{x) = cox~^e~''^ ^ for large enough x had indeterminate 
moments in the above sense. On the other hand, if p{x) were to oscillate around Pfit{x), but sometimes taking much 
smaller values than P^t, then the logarithm will take large negative values; such behavior could lead the integral to 
diverge and allow the moment problem to be determinate. 

To illustrate how delicate the uniqueness issue can be, we note that the probability distribution P{x) — g0(x)e~^ ^ , 
has moments a„ = i(3n + 2)!, that are indeterminate in the Stieltjes sense on [0, oo) by Krein's theorem. However, 
mild modifications of these moments yield determinate problems. For example, by Cor. 4.21 in Ref. 19j, there exists 
a constant c so that the set of moments oq = 1, a„ = c(3n — 1)! is a determinate problem, corresponding to a purely 
discrete probability distribution. 

Overall, we are not able to resolve the question of determinacy, although on balance our expectation is that the 
problem is indeed indeterminate. Certainly we have not been able to find any positive evidence to suggest that the 
moments are determinate. Nonetheless, certain features of the probability distribution can be ascertained. We have 
shown in Appendix A that our moments grow as a power ti mes (3n — 4)!. This rate of growth seems to be just what is 

needed to produce distributions with tails falling as in Eq. (77 1, that is, proportional to a;~^e~°^ ' . We have argued 



that any probability distribution arising from our moments will have a broadly similar tail. This asymptotic behavior 
is all that is needed for many applications of our distributions, such as those discussed in Sect. |VII[ 

It is also worth noting that the conclusion that the probability distribution has a lower bound is independent of any 
concerns about uniqueness, because this follows from existing quantum inequality bounds. Our actual estimates of 
the lower bounds, given in Sect. |IV[ are not strictly independent of the uniqueness issue, but only use a finite number 
of the moments. Thus the numerical answers obtained only depend upon the values of these moments. 



B. Summary 

In this paper we have explored possible probability distributions for averaged quadratic operators in the four- 
dimensional Minkowski vacuum state. We use averaging with a Lorentzian function of time, and investigate the 
distributions for (^^, where (p is a massless scalar field, for ps, the associated scalar field energy density, for E'^, the 
squared electric field, and for pem- In £^11 cases, we infer that the distributions have some features in common with our 
previous results [TS] for a conformal field in two dimensions and for if^ in four dimensions. Specifically, there is a lower 
bound on the distribution, which coincides with the optimal quantum inequality bound on the associated expectation 
value in an arbitrary quantum state. Furthermore, there is no upper bound on the distributions, so arbitrarily large 
positive quantum fiuctuations are possible. 

We have outlined a procedure that, in principle, allows the calculation of an arbitrary number of moments of a 
given distribution. In practice, this procedure can be carried at least as far as the 65*'' moment, which is sufficient to 
allow reasonable numerical estimates of both the lower bounds, and of the asymptotic tail for large argument. These 
are not guaranteed to be unique, but as was argued in the previous subsection, they may be useful. 

If we accept the forms of the tails which we find, then several physically interesting applications follow, including 
the nucleation rates for black holes and Boltzmann brains. It should also be possible to apply these results to the 
study of the small scale structure of four dimensional spacetime, along the lines studied in two dimensions in Ref. j20) . 
It may also be possible to learn more about the non-Gaussian density and gravity wave perturbations in inflationary 
cosmology, which were studied in Refs. |10H12j . Another implication of our form for the tail is that vacuum fluctuations 
will dominate thermal fiuctuations at high energies. The Boltzmann distribution falls exponentially with energy, but 
vacuum energy density fluctuations fall more slowly and hence eventually dominate. 

There is clearly room for further work on the topic of this paper. One obvious problem is to determine whether 
or not the moment problems we have studied are determinate: if so, one would like to know the detailed form of the 
corresponding probability distributions; if not, one would like to know how much information may be extracted from 
the moments, nonetheless, along the lines of the arguments in Sect. [VTl In addition, our results have now trapped the 
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sharp quantum inequality bounds for various operators between the fower bounds given by the methods of Ref. |27| 
and the bounds obtained in Sect. IV which are an order of magnitude smaher. If the moment problem is determinate, 
the latter bounds will coincide with the sharp bound; otherwise it would be interesting to determine what the sharp 
bound actually is. Recall that here we deal only with Lorentzian sampling and only in the time direction. It will also 
be of interest to investigate more general sampling functions, and the effects of sampling in space as well as time. 
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Appendix A: Computation of the moments 



We describe how the moments of smeared Wick squares may be computed for a general derivative 4> of the massless 
field (p in four dimensions, writing p for one more than twice the number of derivatives, so p — \ for '.ip^: and p — "i 
for :ip^:. Thus the two-point function for </>, restricted to the time axis, is given by 

mm'))^^]^ d^^^e— (Ai) 

With smearing along the time axis against smearing function /, the rules for computing the contribution to the n'th 
connected moment of a given connected graph on n vertices may be stated in Fourier space as follows. For each line, 
the form of the two-point function entails that there is a momentum integral over the positive half-line and a factor 
of the p'th power of the momentum; for each vertex there is a factor of f{ujj + ujk) if the vertex is the source of the 
lines carrying momenta ujj and Wfe, a factor of /{cuj — ujk) if the vertex is the source (resp., target) of the line carrying 
momentum ujj (resp., w^), or a factor of f{~ujj — Uf.) if the vertex is the target of the lines carrying momenta ujj and 
ojk] there is an overall factor of (47r^)~" and a combinatorial factor that is 2" for n > 3 and 2 for n — 2. Here / is 
the Fourier transform, defined with the convention 

/oo 
dt/(t)e''^* . (A2) 
-oo 

An important point is that if, as for the Lorentzian, / is real and positive, then every graph contributes positively 
to the moment. Thus any individual graph on n-vertices gives a lower bound on the n'th connected moment. If one 
wishes to compute the dimensionless moments, defined in the text so that a„ = (47rT(P+^)/^)^"/i„, the overall factor 
2"(47r2)-" is replaced by 8'V-(p+i) for n > 3 (or 32 in the case n = 2). 

In the particular case of the Lorentzian (|9]), we have f{Lu) = e"'"'"^, and a simplification of the computation rules: 
a vertex met by lines carrying momenta LOj and ojfc contributes Q~\^i~'^''\r ^ target for one and a source for the 

other, or ^-^i'^i+'^k)T otherwise. This means that the overall integral over wi, . . . , w„ factorizes at each vertex that is 
either a double source or a double target. 

Recall that the graphs involved are drawn on n vertices xi, . . . , a;„, placed in increasing order from left to right. 
Each vertex is met by two distinct lines, and each line is directed to the right. In particular, xi is the source of both 
lines connected to it. We may represent such a graph by a permutation a of the set {1, . . . , n} of integers, subject to 
the conditions that cr(l) — 1 and cr(2) < ct(7i). To reconstruct a graph from a permutation, draw lines from xi — x^^i-^ 
to Xcr(2), from x„(2) to a;cr(3); and so on, finishing with an line from x^r^n) to xi. [For this reason, it is convenient to 
adopt a convention that cr(n -|- 1) = 1.] Then place rightwards-pointing arrows on each line. On the other hand, to 
encode a graph as a permutation, start at xi and follow the shorter of the two lines to the vertex it meets [i.e., of 
the two vertices joined to xi, choose the one with the smaller label], and continue along the other line meeting that 
vertex. At each subsequent vertex, continue along the line not previously traversed, eventually returning to xi. Then 
(j{k) is defined to be the fc'th vertex met on this trip. 

A run of cr is a set of consecutive integers in {1, . . . , n -|- 1}, say p,p -I- 1, . . . 9, such that cr(p), a{p -I- 1), ... , (j{q) is 
a monotone sequence, either ascending or descending, and so that no superset of consecutive integers has the same 
property. The length of the run is defined to be g — p. Every permutation used to label our graphs corresponds to an 
even number of runs, alternating between ascending and descending, whose lengths sum to n, and with consecutive 
runs sharing a common endpoint. For example, the permutation 14536782 (i.e., (t(1) = 1, cr(2) = 4, (7(3) = 5,...) has 
runs 1,2,3; 3,4; 4,5,6,7 and 7,8,9, of lengths 2,1,3,2; representing each run by its image under the permutation, 
these runs are more transparently written as 145, 53, 3678, 821. 
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The contribution to the n'th connected moment arising from the graph corresponding to any given permutation is 
easily seen to factorize into terms corresponding to the runs, whose values depend only on the length of the run: in 
our example, the graph contributes ^^K2KiK^K2 = %^KiK2K'i to the dimensionless connected moment Cg, where 
the Kj correspond to the special case Kj — Kj'^'' of the family of integrals 

Ki['^ = ^ / dkidk2--- dfc„/cf+'^(/c2 • • • fc„)Pe~'=ie-^-i' Ifc.+i-fcdg-fc". (A3) 

[Here, ki = ujiT are dimensionless versions of the momenta previously used.] 

These considerations reduce the computation of the n'th connected moment to two problems: the computation of 
the Kj and the enumeration of all permutations in the class considered with a given run structure. To address the 
first of these, we note the easily proved identity 



dfcFe-e-l-'^l = |^$:i^, (A4) 

=0 



of which the standard integral 



2P+ 

is the K = special case, and which entails the recurrence relation 



*''e-» = -5pj (A5) 



p+r+1 



r'=0 

As Af'^ = 2-(P+i)p!(P+''), it follows that K^P is given by 

^ / \ / r„_i=0 r„_2=0 1-1=0 A;=l \ ^ ^ 

for any integers n > 1 and r > 0. Although we have not found a closed form expression for the Kn \ the above 
expressions allow for them to be computed efficiently. 

To the best of our knowledge, the problem of enumerating permutations of the class we study in terms of their 
run structure does not appear to have been solved in the literature on enumerative combinatorics, although related 
problems have been studied for over a century. A closed form answer seems out of reach, but generating function 
techniques allow one to build up a solution for each n in a recursive way. The details will be reported elsewhere 137) , 
but the overall result is the following: for each n, let /C„ be a polynomial in the variables Ki, . . . , i^n-i, with JC2 = 
and subject to the recurrence relation 



dlCn-1 . dJCn- 



Then, for n > 3, the coefficient of K^^ ■ ■ ■ K^^^^ in /C„ is precisely the number of permutations a of {1, . . . , rt} with 
mi runs of length ^ (1 < ^ < n — 1), subject to the side conditions cr(l) — 1, cr(2) < a{n). In the case rt = 2, we find 
half of the number of such permutations. 

The generating function is extremely convenient, because it already incorporates the sum over all possible connected 
graphs. Putting this together with the other considerations above, the n'th dimensionless connected moment is given 
by C„ = 8"/C„, for any n > 2, where the variables Kj are given the values defined above by (A3) (recalling that 
Kj = K^j^'^). For example, we find the explicit formulae: 



C2 = ?>2Kl (A9) 

C3 = S^i^ai^i (AlO) 

C4 = 8^ (i^3^i + Kl + K*) (All) 

C5 = 8^ {K^Ki + 3K3K2 + SK2Kf) (A12) 

Ce = 8^ {K5K1 + ml + 'iKiK2 + liK^Kl + HKIkI + 8Kf) (A13) 

C7 = 8^ (KeKi + IOK4K3 + 5K5K2 + 19KiKf + eGK^Ki + 123K3KfK2 + IS6K2KI) (A14) 
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which can be used to provide the first few connected moments for -.Lp^: in the case p = 1 or :(^^: in the case p = 3. 
One may check that the coefficients inside each parenthesis sum to [n — l)!/2, the total number of connected graphs 
involved in the n'th moment. 



Appendix B: Asymptotics of the moments 

In this appendix we give asymptotic estimates for the n'th moments of the Lorentzian smearing of the Wick square 
of the \{p — l)'th derivative oi if as n becomes large. We rigorously establish a lower bound and also give an upper 
bound, for which our reasoning is not completely rigorous, but which appears to be satisfied on the grounds of 
numerical evidence. The basic observation is that the dominant contribution to C„ [and hence the full dimensionless 
moment o„] is this is certainly a lower bound (as all terms are positive) and numerical evidence suggests 

that it gives the correct answer modulo a fractional error of order n~^. Thus lower bounds on the Kj will give rigorous 
lower bounds on C„, while upper bounds give an upper bound on the C„ that seems reasonably secure, albeit not 
fully rigorous. In terms of permutations and graphs, the dominant contribution arises from the identity permutation, 
and thus the graph on n vertices that has lines from Xk to x^+i for each fc = 1, . . . , n — 1 and an line from xi to x„. 
The graphs in Fig. [l] represent the case n = 2 and n — i. 

We begin with the lower bound, which is 

K^:'>{''^'^'}:'V^i^^\ (Bi) 



where 



n{p+l) - 1 



^ (n(p+l))! \^ r + n[p+l) 

n!(2P+i(p+l))" r + fc(p+l)' ^ ' 



in which the product over k is taken to be equal to unity in the case n — 1. The bound (Bl I is proved by induction. 



noting that the it is true (indeed, an equality) for n = 1. Supposing that it holds for some n > 1, we use the recurrence 



relation Eq. ( A6 ) to show that 



r^ir) > ^ (p + r\({n + l){p+l) + r\ (^+1+,) 
"+i-2P+ilv P )\ n{p+l) j^" 

_ p\ r + (n + l)(p+ 1) /(n + l)(p+ 1) - 1\ /(?i + l)(p+ 1) - 1 + 



2P+1 r+p+1 V p J\ (n+l)(p+l)-l 



LiP+^+r)^ (B3) 



(r) 

where, in the first line, we have used the fact that the constants L„ are clearly monotone decreasing in r for each 
fixed n, and the identity (0.151.1 in Ref. [38 ) 

the second line is an elementary algebraic manipulation. A further algebraic manipulation shows that 

r{r) _ pl r+{n + l){p + 1) /(n + l)(p + 1) - 1\ 

K+i-^ ^T^Tl ^ p j " ^ ' 

(r) 

which allows us to conclude that the bound on Kn holds for all n by induction. Noting that 

r(o)^ (n(p+l))!n" 

(2(P+i)(p+l))»(n!)2 ^ ' 



we obtain a lower bound on J„ = K^^ Kf\ = p!2 ^p^^^'K^^^ as 



J >^r(0) ^ (p+l)!((n-l)(p + l))!(n-l)"-^ 

« - 2P+1 (2P+i(p+l))"((n-l)!)2 ■ ^ ' 
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In a similar way, we find an upper bound 
where 



n{p ' " 



" n! (2P+1 (p +1))" ^J-^ ^J-^kp + r + n + q- 



1 



(B8) 



(B9) 



and the product on k is again regarded as a factor of unity in the case n = 1. From this expression, it is clear that the 

(r) 

are monotone increasing in r for each fixed n. The double product can be also be written as a ratio of products 

case r = tt 

k{p+l) + q 



of F-functions and other simple functions; in the case r — Q there is a particularly simple expression 

n—2 p 



nn 

k=0 9=1 



kp + n + q — I 



(BIO) 



As before, we prove (B8) by induction, noting that it holds with equality in the case n — 1. Supposing that it is true 
for some n > 1, the recurrence relation Eq. (A6| gives 



r^ir) ^! fp + r 



fp + r\ "t^' /n(p + 2)-2 + r'\ 
'\ P J .^A n{p+l)-l ) " ■ 



(Bll) 



Over the summation range, we have Un •* < u'^^^^^\ so 



r'=0 



n(p + 2) - 2 + r' 
n{p + 1) - 1 



l+r) 



. , p+r+n y , . 

p + r\ pi \^ f n(p + 1) — 1 + r 



p J 2P+ 



r E 

r"=0 



n{p+l) - 1 



(B12) 

where we have changed summation variable to r" = r' + n — 1 and extended the summation range in the second step. 
Using Eq. ( B4 1 , this gives 



M ^ _p}_fp + r\ (n{p + 2) + P + r\^(p+i+r) 



^ 2P+1 



n{p+l) 



Using the fact that 



p + r\ / n{p + 2) + p + r 
n{p + 1) 



{{n + l){p + /(n + l)(p + 2) -2 + r\ -A- r + q 



[n + l){p + l)[n{p + l)]\p\ \ + 



^ z; - z ^ ; \ -pr 



r + n -\- q 



(B13) 



(B14) 



it is then easy to show that ( B8 ) holds with n replaced by n + 1 and hence for all n by induction. 
We may then obtain the upper bound on ,/„ as 



Jn < 



(p+l)!(p+l)3((n-l)(p + 2)-2)! 
(2P+i(p+l)2)" ((n-2)!)2 



(B15) 



after some manipulation. 

Using Stirling's formula, (nA — B)l ^ ^/2n{nA/e)'^^~^~^^^^e'~-^^^^'^ . Then one may check that the lower bound 
in (B7| is, asymptotically. 



Jn > 



{p + iy. f p 



p+1/2 



27re \p + 1 

while a similar calculation at the upper bound gives 



{p + Ife 
pP2P+i 



(np- (p+1))! 



P 



{p + iy.{p + iy 

27r(p + 2)3 \p + 2 



p+l/2 



[p + 2)P+2 
2P+i(p+ 1)V 



{np-{p+l))\ 



(B16) 



(B17) 
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so the ratio of the upper bound to the lower bound grows as ~ a/?" as n — cx3, with 



P + 2J e \p+l 

which in the case p = 3 gives a = 0.6373520649, /3 = 1.122678959. So we have a reasonable control over the leading 
order contribution. 

As mentioned above, it is certain that the dimensionless moment a„ obeys a„ > 8" J„, and numerical evidence 
suggests that a„ ~ 8" J„ at least in the case p — 3 (we believe that this is true for all p > 1 and could be proved with 
more effort). On that basis, we have 

0.513395 X 3.221667" < ^g^"" < 0.327213 x 3.616898" (B18) 
in the p = 3 case, for n oo. This supports the growth estimates given in the text. 
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